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Abstract 


The recent fusion of decades of advancements in mathematical models, numerical 
algorithms and curve fitting techniques marked the beginning of a new era in the science 
of simulation. It is becoming indispensable to the study of rockets and aerospace 
analysis. In pneumatic system, which is the main focus of this paper, particular emphasis 
will be placed on the efforts of compressible flow in Attitude Control System of sounding 
rocket. 


Key words: Sounding rocket, mathematical model of pneumatic system, pressure PDE, Euler finite 

difference . Cubic Bezier, MatLab 


Introduction 


The development of an Attitude Control System (ACS) is very complex due to the 
nonlinear nature of the expressions of the system. Pneumatic analysis is even more complicated 
by the fact that the operating fluid is extremely compressible. Most problems which occur with 
pneumatic applications become apparent in the dynamic operation. To analyze the system during 
dynamic operation it is necessary to develop component and system models. This paper 
discusses the modeling approach and reveals the principles upon which it is based. The technique 
unifies the various disciplines to effectively design and also modify existing design. It formulates 
the governing dynamic equations based upon the topological information presented by the 
schematic. The integrity of the output is cross validated by comparing to empirical data and real 
gas properties extracted from NIST (National Institute of Standards and Technology). 


Pneumatic Schematic of Sounding Rocket 

The pneumatic schematic consists of a high pressure reservoir, a manifold, pitch and roll 
blocks. It is divided into three segments, the rationale that schematic is divided into three is 
because of the provided experimental data at the end of each segment. 

• Segment 1 (form tank to regulator) 

• Segment 2 (from regulator transducer to pitch block) 

• Segment 3 (from regulator transducer to roll block) 


Mathematical Model of the Pneumatic System 

Introduction 

In this analysis, a descriptive model of the system is built as a hypothesis of how the 
system could operate, and to estimate how the design variables could affect the system. 

The mathematical model describes a flow system by a set of fluid dynamics variables 
and a set of governing equations. Applied equations establish relationships between the 
flow variables. The variables represent fluid properties of the system. The equations are 
derived based on fundamental physical principles. 


l 



Roadmap for Deriving Nonlinear ODE for Pneumatic System of Sounding Rocket: 



Fundamental physical 


principles 


Real model of 
the flow 





Nonlinear 
Differential 
Equations 
Ordinary in time 


Governing equations 
of fluid flow 


Bernoulli’s equation 


2k sk 

Fluid Dynamics 
Model of the system 


Continuity equation 
Real gas 

Van Der Waals eq. 




Assumptions and 
concepts 



Isentropic 
Adiabatic 
No shock wave 


Steady Flow 
Homogeneous 


Fig 4.1 


Modeling technique 

Mathematical model of pneumatic system is subcategorized into blocks with 
inputs and outputs. Each block represents a pneumatic component in the assembly. A 
posteriori pressure behavior is provided from experimental data. Since there is no priori 
data available for critical blocks, the posteriori data is utilized to simulate the behavior. 
Therefore the model consists of a set of inputs, applied governing equations and a set 
of outputs. 


Gas flow through the system (assumptions and concepts) 

The analysis of gas flow in the system involves a number of assumption and concepts: 

• The density of the fluid changes with respect to pressure v flow is considered to 
be a compressible flow. 
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• The gas expands isentropically (i.e., at constant entropy) v the flow is reversible 
(frictionless and no dissipative losses), and adiabatic (changes in temperature 
occur due to changes in pressure of a gas while not adding or subtracting any 
heat.) 

• The gas is assumed to be a real gas (Van Der Waals equation) and 
homogeneous. 

• The system of flow through the tubes could be in a steady state because there is 
a constant flow of fluid. Conversely, the tank which is being drained or filled with 
fluid would be a system in transient state, because the volume of fluid contained 
in it changes with time. 

Symbols 


Symbol 

Definition 

V 

Tank Volume 

u t 

Velocity of flow gets in 

u r 

Velocity of flow when gets out (downstream) 

Ui 

Velocity of flow in component 

A, 

Area flow gets in 

A r 

Area flow gets out 

Li 

Length of the tube 

Di 

Diameter of the tube 

Z 

Elevation 

€ 

Surface roughness of the tube 

Ki 

Loss coefficient 

fi 

Friction loss 

P 

Density of the fluid inside tank 

771 

Mass flow rate 

G 

Gravity 

R s 

Gas constant 

T 

Absolute temperature 

Pt 

Tank Pressure 

Pr 

Regulator Pressure 

N 

Number of moles 

A 

Measure of the attraction between particles 

B 

The volume excluded by a mole of particles 

Y 

Specific heat dimensionless ratio 

Y, 

Expansion factor 


Table 3.1 
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Governing equations: 

In this section: 

• The fundamental physical principles are written down 

• Then applied to the suitable model of the flow. 

• Obtain equations which represent the behavior of the flow. 


1. Continuity equation 

Consider the flow model i.e. (segmentl from tank to regulator), namely, a control 
volume of arbitrary shape and a finite size. Conservation of mass states: 

“Net mass flow out of control volume through surface is equal the time 
rate of decrease of mass inside control volume.” 

dp 

£ + v.(p^) = o 


Assume the gas is homogeneous v 


Given an area A, and an incompressible fluid flowing through it with uniform velocity U 
with an angle 0 away from the perpendicular direction to A, the flow rate is: 

Q = A. U. cos 9 

In the special case where the flow is perpendicular to the area A, that is, 0 = 0, the 
volumetric flow rate is: 

Q = AU = A t U t = A r U r ^U t = (^) U r eq. (a) 

Flow expands isentropically, the flow rate is derived for incompressible flow then it is 
modified by introducing the expansion factor Y t to account for the compressibility of 
gases, v 


m = Y i .p.V = Y i .p.Q eq.(b) 
From eq.(a) and eq.(b) v nv = Yi.p.A r . U r eq. (c) 
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The expansion factor Y t , which allows for the change in the density of an ideal gas as it 
expands isentropically is given by: (Ref: Perry's Chemical Engineers' Handbook 7 th 
edition) 


Y t = 



2. Equation of state (Real Gas) 

The Van Der Waals equation of state for a fluid such as Argon is: 

( P + ^)-(V-nb) = n.R,T 

By dividing the side of equation by the mass of the gas, the volume becomes the 
specific volume. 


f n 2 a\ /V nb\ n. R. T 
\ + V 2 ”) ' (in “ "m ) = 


m 


t 


n 2 a\ /V nb\ 

P + -T75- . = R S .T 

V 2 J Vm m / 

( n 2 a\ / nb\ 

(P+ ) = R s- T 

/ n 2 a\ / , nb\ 

( p + V 5 ") ' ( P _ nr) = R *' T 




nb 

m 


( n 2 a 

( P+ VT 


2 a> _1 


,R S .T 


nb 

H 

m 


p-'=( p +^) - r s - t 

P -(^r, s ,^y e q . ( d) 


5 



Absolute temperature refers to use of the Rankine (°R) temperature scales, with zero 
being absolute zero. The value of new constant depends on the type of gas as opposed 
to the universal constant gas which is same for all gases. 


dp 

dt 



dp_ 

dt 


d_ 

dt 







eq.(e) 


3. Bernoulli’s equation 

In the present section, the third physical principle as itemized in the roadmap applied to 
the pneumatic system, namely, 

“Energy is conserved; meaning rate of change of energy inside fluid element is equal to 
the total of net flux of heat into element and rate of work done on element due to body 

and surface forces. ” 


In general form the energy equation is derived from Green Gauss’ theorem as following: 

U 2 \ ff P U 2 . . 

u + — + gz) pdV + jj> p(u + — + — + 9 Z )V- = Q - W s 


By assuming steady state for the pneumatic system: 

P U 2 - . . 

p(u H I- — + gz)V. nd.A = Q -W s 

P 2 




Furthermore for compressible, adiabatic flow of perfect gas ( Pp Y = constant), steady 
flow and negligible shock waves, the Bernoulli’s equation represents the simplified 
energy equation as following: 


y-i' 


U r 2 -U t 2 y P t l . fP r \Y 


2 y 


Y ~lj 1 "(^) r ) + 9 (z r -z t ) = 0 


Y- IN 


u r 2 -u 2 Y Ptl. fPr\y 
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By substituting velocity equation U r and density p from real gas law (eq. d) into derived 
mass flow rate (eq. c) ••• 


m r = Yi.p.A r . U r 



eq.(f) 


By assuming the tank as the control volume and the gas as homogenous object v 

rn t = Y i j t (pV) = Y t V ^ eq. (g) 

By substituting changes of density in time eq.(e) into eq.(g) v 


= r t £(pV) = Y i V^ = Y l .V.R s ~ l .f t [T^- 1 (p w + ^)] eq. (h) 


Principle of mass/matter conservation states that the mass of a closed system (in 
the sense of a completely isolated system) will remain constant over time, v 


m t = riir = m eq.(i) 

Therefore mass flow rates calculated from eq.(h) and eq(f) have to be equal, v 
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d f . / n 2 a\ 

di [ T(t) \ P V + ^P) 


= R s .p.V~ 1 .A r . 

y 

-[( 
p ' 


- 2 g(z r - z t ) 


Eq.(j) (Partial differential equation for first segment of the pneumatic system) 
Classification of derived differential equation 

The mathematical model of the pneumatic system is classified as of the following: 

1. Linear vs. Nonlinear: Since the operators in this mathematical model exhibit 
nonlinearity, the resulting mathematical model is defined as nonlinear. Although linearity 
exists in the system. 

2. Deterministic vs. probabilistic (stochastic): The applied model is one in which every 
set of variable states is uniquely determined by parameters in the model and by sets of 
previous states of these variables, therefore it is deterministic. 

3. Static vs. dynamic: The curve of pressure drop is descretized in time, and it is dynamic 
since a dynamic model accounts for the element of time. Because it is a dynamic model, 
it is represented with differential equations. 

4. Lumped vs. distributed parameters: Argon (Ar) is the operational fluid in this model. It 
is homogeneous (consistent state throughout the entire system), then the parameter is 
classified as lumped. 


Empirical Data of the Pneumatic System of Sounding Rocket 

Pressure test: 

Transducers and gages are utilized to measure the pressure in the test article, 
for procedure please reference to appendix B. 

Thermal Testing Procedure: 

Thermal testing for the system began with calibrating each of the thermal 
couples. This was achieved using a calibrated thermal chamber along with two 
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calibrated thermal couples. The calibration was done using three temperatures, 20 
degrees C, -60 degrees C and 40 degrees C to achieve a calibration curve. 


Once the calibration curve was known testing could begin on the pneumatic 
system. The thermal couples were placed in or around the flow of the system to 
determine the temperature of the flow. After placing the thermal couples in the desired 
location the system was pressurized. Once the system was stable the thermal data 
acquisition system was turned on and the pitch valves of the system were opened until 
the pressure transducer in the system read 3700 psi. The test was complete once the 
pressure stabilized at 3700 psi for nine-hundred seconds. 


Temperature of Argon in Tank from Experiment 



Fig 5.1 

Computational Fluid Dynamics Approach 

Introduction 

In this section, computational fluid dynamics approach synergistically complements the 
other two approaches of pure theory and pure experiment. 
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The system of differential equations for n segments that represents the pneumatic 
system: 


r/UOi 


fl(t,P t ,P r , ... 

-.Pn) 1 

p r (t) 

— 

f 2 (t,P t ,P r , - 

■ ’iP-n) 



.f n (t, P t , P r , ... 

• • i P-n) - 


The Jacobian is an n-by-n matrix of partial derivatives: 


\dfi 

dfi 

dfi ] 

dP t 

dP r 

.. dP n 

df 2 

df 2 

. df 2 

dp t 

dP r 

; dPn 

dfn 

dfn 

■ df n 

ldP t 

dP r 

dp„\ 
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The influence of the Jacobian on the local behavior is determined by the solution to the 
linear system of ordinary differential equations generated for each segment of 
pneumatic system (Eigenvalue problem): 


Derived pressure differential equation (eq. j) is a first order nonlinear partial differential 
equation. By application of temperature test results, and injective concept (If Pressure 
and Temperature both are injective, then P °T is injective.) Therefore it converts to 
nonlinear ODE in time with initial value problem: 


^ = /(t,P t (t),P r (t)) f:R*R -» fl 

P ttlme.0’ P rtime-0 = Siven 

Because of the non-linearity, it cannot be solved analytically and exactly. Since the 
regulator pressure is discretized in time domain, a MATLAB iterated program is 
developed for the solution. The solution is an approximation, based on forward 
difference explicit Euler method. It uses a step size and generates the approximate 
solution. The smaller the time step the more accurate the solution. 

P t (t) = ^(r+ft) p t(r) -> p f (t + h) = P t (t ) + hf(t, P t ) ( numerical solution ) 

t(n + 1) = t(n) + h\ 


However, the exact solution shall be calculated by Taylor series. Taylor series of 
P t (t + h ) that is infinitely differentiable in a neighborhood of a real (t) or complex (t), is 
the power series: 

P t (t + h)= P t (t ) + hf(t, P t ) + < tfc+i ( exact solution ) 

But P c (Cfc) is not known (no exact solution is available), by comparing the exact solution 

h 2 

and numerical solution, the —P t (£ fe ) is truncation error. It is second order of the 

associated time step 0(h 2 ). The other error is round-off calculation error arising from the 
use of floating point arithmetic. 
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The Bezier Curve Smoothing Technique: 

Extracted upstream pressure data in time domain encountered truncation error and 
round off error. Based on iterated expectation law, we need to smooth the fitted curve. 
Parametric Bezier curves are important tools used to model smooth curves. In this study 
Cubic Bezier fitting method is applied to generated upstream pressure in attitude control 
system. Cubic Bezier has four control points and the parametric form of the curve is: 

B(t) = (1 - t) n P 0 + Q (1 - t) n ~ 1 tP 1 + Q (1 - t) n ~ 2 t 2 P 2 + - + Q t n P n ; t e [0,1] 

n = 3 ( Cubic Bezier ) -> B(t ) = (1 — t) 3 P 0 + 3(1 — t) 2 tP 1 + 3(1 — t)t 2 P 2 + t 3 P 3 ; 
t E [0,1] P 0 , P lt P 2 , P 3 are Bezier Control Points. 

Flow through Segment 1 (from tank to regulator) 

In the pneumatic system, flow of gas from tank to the regulator is considered as 
segment 1 . The objective is the calculation of pressure, mass flow rate and velocity at 
downstream. This segment consists of a tank, a high pressure block with expansion and 
contraction, a regulator and an adiabatic flow through constant diameter tubes with 
friction and bends, (assume volume, temperature and pressure of the tank, regulating 
pressure, tube bend angles, lengths diameter and surface roughness are known from 
experimental data.) 

Procedure: 

1. Calculation of friction factor for the tubes: 

• Estimate Velocity 

• Compute Re = ^ and then / from Colebrook-White equation 


/ = 


1.14 - 2 log(- + 


e 21.25 


-2 


D Re 


.9 


• Calculate U from U = (^-^) / 

V fpL ) 

• Iterate, returning to step 2, until desired accuracy is achieved. 

The coding in MatLab: 

for U_est=.l: .25:1000 
Re ynold_comp=U_e s t* D_comp / Ar_neu ; 

f rictionl= (1 . 14-2* loglO (epcilon/D_comp+21 .25/ (Reynold_comp A . 9) ) ) A -2; 
de 1 t a _P _c omp 1 = P_c omp 1 - P_c omp 2 ; 
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U_comp= (2* delta_P_compl* D_comp/ (f rictionl* Rou_var_Ar* L_comp) ) A . 5; 

if (U_comp <= (U_est + err_rate) & U_comp >= (U_est -err_rate) ) 

break 

end 


end 


2 . 

3. 

4. 


5. 


Input loss coefficients for tube bends from Fluid Dynamics handbook. 
Calculate losses in abrupt contraction and expansions from Applied Fluid 
Dynamics handbook. 

Calculate velocity at upstream from: 


u t = 




Z n fi^i I ym jz _ Y 1 

i= l D t 2 y 




Solve the nonlinear partial pressure drop differential equation: 


d_ 

dt 




r s . P y-\A r 


2 P t 

u t 2 + — 

p 



- 1g(z r - z t ) 


Pseudocode in Matlab for Euler method: 


function 

[ Matrix] =Euler_ODE_solver (n, d_orifice f Area_Orifice, Volume_Tank f gr f Temp f R , A , Di 
scharge_coef f , Rou_var, z_2 f z_l , U_t , gamma) 
for i=2 : length (A) 

if (Tank_Pressure ( i — 1 ) -Tank_Pressure_dif ferential (i-1 ) ) > A(i f 2) 

Tank_Pressure (i) =Tank_Pressure (i-1) -Tank_Pressure_dif ferential (i-1) ; 

else 

Tank_Pressure (i) =A (i f 2) ; 

end 

if complex_root_checker (i) <0 | shock_critical==Ratio_Mach 
break 
end 

Tank_Pressure_dif ferential (i) =Tank_Pressure_dot (i) .* (A (i f 1) -A ( (i-1) , 1) ) ; 
Delta_P (i) =Tank_Pressure (i) -A (i f 2 ) ; 

end 

[ Bezier_generated_data] =NGC_All_Bezier_Interpolated_values (ctrlPO f ctrlPl , Ctrl 
P2 f ctrlP3 ) 

[ ctrlPO r ctrlPl , ctrlP2 r ctrlP3 f fbi] =NGC_bezier (DeltaP_f rom_venting f Max_Square_d 
istance) ; 

[ Vel,mdot] =mdot_calc (Pressure_vs_Time f M ( : f 2 ) , Rou_var_Ar f Area_Orifice f U_t f gamm 
a_Ar f z_l , z_2 , gr ) 

6. Calculated mass flow rate is constant for the system, it is used for segment 2 and 
3 to calculate the pressure drop at roll and pitch blocks. 
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Flowchart of the set of the MatLab programs (Seg_1) 


Reads 

Reg. 

Pressure 


ColeBrook White friction 
loss calculator 


^EstimatecP^Yes 
Velocity is — • 

Xaccurate/^ 


Subsonic Compressible 
Velocity, tubejoss 
Bendjoss 
calculator 


Euler Nonlinear Solver 


/ Anomaly^ 
chkr: shock, 
^x>mplex rooj> 


Output: 

Regulator 

Pressure 


Bezier Fitting and Smoothing Technique 


Modify 

No > 

/Max 

Euclidian distance 

Bezier control 

inputs 

* 

. distance 

index for segments 

points generator 


Bezier 

Interpolated 

values 


Cubic 
Bezier Ctrl 
Points 


All Bezier smoothed 
interpolated values 


Output: 
Smoothed 
Regulator P 


X X — 

Empirical vs. Theoretical 

Yes 

^ / Within \ p- 

Output generator 

comparator 

^Njange^/ 

call 

| Modify [ i 

No 

v 


Fluid Dynamics 
data 

Vel,mdot,P 


Plot Pressure, 
mdot vs time 





Flow through segment 2 (from regulator to roll block) 

In the pneumatic system, subsonic flow of gas from manifold to the roll block is 
considered as segment 2. This segment consists of an ideal adiabatic expansion or 
contraction at changes in manifold and adiabatic flow through constant diameter tubes 
with friction). Based on calculated mass flow rate and fluid dynamics properties at 
regulator (upstream) from segment 1 , now the pressure, density, velocity and 
temperature can be calculated at roll block (downstream). 

Procedure: 


1. 

2 . 

3. 

4. 


5. 


6 . 


7. 


The mass rate of flow through each section is constant; 

The isentropic stagnation temperature and isentropic stagnation speed of sound 
are constant for adiabatic flow regardless of friction. 

The isentropic stagnation temperature, pressure, density and speed of sound are 
constant across isentropic expansion and these values varies only with Mach 
number and the ratio of specific heats. 

Using conservation of mass: 


rhR 1 / 2 T 0 1/2 / 

~~P^ = Y M { 1 + 


7-1 


M 2 


y+1 
*2(i -y) 


The ratios of static pressure, temperature, and speed of sound to their isentropic 

stagnation values are functions of Mach number and ratio of specific heats. 

T c. „ ( y ~ 1 9 

1 +-^^M 2 


= (— ) 2 = ll+^—^M 2 ) 

T 0 c 0 \ 2 / 


f ( y-i ,y- 

_ = ( 1+ _ m2 ) 


Based on solution to derived equation at regulator (segment 1 ), the followings 
are known: 

The mass flow rate 
Static pressure 
Absolute temperature 
Density 
Velocity 

Speed of sound (from NIST table) 

Ratio of specific heats(from NIST table) 

This implies that M = ^ and value of ^ is derived from following formula for 


adiabatic flow: 


fL u 1-M 2 y — 1 (7 + l)M 2 

D yM 2 2y ln 2[l + (y - l)M 2 /2] 
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8. Friction factor for the tube is calculated by Colebrook-White equation and as a 
function of surface roughness and Reynolds number: 


/ = 


6 21.25 


—2 


114 - 21o % + ^ 


9. The frictional losses between two points are the sum of the straight pipe loss, 
bends loss, dividing T loss, valve loss, expansion loss and manifold loss. 

(fL 12 


(l— 11 
V Di 


T K bends "F ^mnfld K v +Kf 


10. Then from 


f L M 7 


Dt D 1 

one can calculate Mach number by interpolation method: 

Mach2] =Mach_interp_calc (fL_M2overD, gamma) 


D i 

function 


Friction Length vs Mach for Adiabatic Flow of Argon Through Tubes 



Friction Length vs Mach for Adiabatic Flow of Nitrogen Through Tubes 



Fig 9.1 

T c 

11. The ratios of — and — can be calculated from: 

T o c 0 
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T 

% 



y — i 

2 



) 


-i 


12. Since the isentropic stagnation temperature and speed of sound are constant 
across the pipe: 

T2 = T °( T / t 0 ) 

C 2 = C o( C /c 0 ) 

13. From M 2 the flow velocity at 2 is known now U 2 = M 2 c 2 . 

14. Conservation of mass implies that the density must be p 2 = 


15. The final step is the calculation of Pressure and Thrust at the end (Van Der 
Waals equation): 

P2 A 
P = p x — -f 
Pi A 


Thrust = j> p.h.dA + j) p. V. (V. n)dA 


Pseudo code in MatLab for roll segment: 


fL_MloverD= (l-Machl A 2) / (gamma* Machl A 2 ) + ( gamma - 

1) / (2* gamma)* log ( ( (gamma+1 ) *Machl A 2 ) / (2* (1+ (gamma-1) *Machl A 2/2) ) ) 
ratTloverT_0= (1+ (gamma-1) /2*Machl A 2) A -1 
ratCloverC_0=sqrt (ratTloverT_0 ) 

T_R_0=ratTloverT_0 A -l* T_r 
c_R_0=ratCloverC_0 A -l* c 

K_losses=K_mnf ld+nT* K_T+nb90* K_90bend+K_roll_preece_valve+K_ctrction 
f L_M2overD=fL_MloverD- (fLoverD+K_losses ) 

[ Mach2] =Mach_interp_calc (fL_M2overD f gamma) 
ratT2overT_0= (1+ (gamma-1) /2*Mach2 A 2) A -1 
T_2=ratT2overT_0* T_R_0 
ratC2overC_0=sqrt (ratT2overT_0 ) 
c2=ratC2overC_0* c_R_0 
U_roll=Mach2* c2 
rou2=mdot_system/ (U_roll*A) ; 

P_roll=P_r* T_2/T_r*rou2/roul 
Thrus t=P* A+rou2* U roll A 2*A 
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Flowchart of the MatLab programs (Seg_2) 



Fig 9.2 


Flow through segment 3 (from regulator to pitch block) 

In the pneumatic system, subsonic flow of gas from manifold to the pitch block is 
considered as segment 3. The procedure for pressure drop analysis is similar to 
segment 2, but the characteristics of the components are different. 

Pseudocode in MatLab for pitch segment: 


roul=mdot_system/ (U_r*A) 

fL_MloverD= (l-Machl A 2) / (gamma* Machl A 2 ) + ( gamma - 

1) / (2* gamma)* log ( ( (gamma+1) *Machl A 2) / (2* (1+ (gamma-1) *Machl A 2/2) ) ) 

[ K_mnfld, K_90bend] =pitch_loss_coef f (qty_mnfd, nb90) 

K_losses=K_mnfld+nb90*K_90bend+K_ctrctionfL_M2overD=fL_MloverD- 
( f Lover D+K_los ses ) 

[ Mach2] =Mach_interp_calc (fL_M2overD, gamma) 

P_pitch=P_r* T_2/T_r* rou2/roul 
Thrus t=P_pitch* A+rou2* U_roll A 2* A 
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Flowchart of the MatLab programs (Seg_3) 



Fig 10.1 

Evaluation of Applied Model 

In this section the followings are addressed: 

• How does the pneumatic analysis program work? 

A. Segment 1 (From tank to regulator) 

1. Inputs: 

i. Characteristics of the components of the system (i.e. Volume, 

Diameter, Length, bend angle, surface roughness, etc.) 



Command Window 



_j 


^ To get started, select MATLAB Help or Demos From the Help menu. 


Regulator Fraction Fully Open : 1 
Enter Volume of the tank in cubic inches: 200 
Enter quantity of 45 degrees bend: 3 

Fig 11.1 Command window in MatLab 
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ii. Operating temperature and regulating pressure 


Pressure_vs_Time=xlsread ( ' C :\ hyekakan\ j armen\ imput . xlsx) 


iii. Properties of fluid for specified pressure and temperature 
(state of the gas) by referencing NIST software. 


REFPROP (argon) - NIST Reference Fluid Properties 


File Edit Options Substance Calculate Plot Window Help Cautions 



Fig 11.2 Generated outputs in NIST 


density of Argon in Tank 



Fig 11.3 

2. MatLab files: 


A set of programs is developed to model the pneumatic system. By 
running main file and above inputs, it is linked to appropriate .m files 
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(function call) and each .m file does its task and returns outputs to the 

main file. 

i. main.m (is where the program reads inputs and starts 
execution) 

ii. Velocityandl_oss_calc.m(iteratively calculates the velocity in 
upstream and losses associated with the system from tank until 
manifold) 

iii. Euler_PDE_solver.m (solves derived partial differential 
pressure drop equation for the pneumatic system, output 
is either deltaP or pressure in upstream) 

iv. NGC_bezier.m (Approximation of data by Cubic Bezier 
curves. Based on least square fit, uniform 
parameterization. Finds Control Point of Bezier Curve that 
approximates the given data up to specified squared 
distance limit 

v. NGC_columnVec_output.m(changes row to column 
vector) 

vi. NGC_Euclidean_distanceandlndex_getter.m (this 
algorithm is based on Euclidean distance) 

vii. NGC_vector_checker.m (checks the argument is a vector) 

viii. NGC_Bezier_AIICtrlPoints_generator.m(this program finds 
the control points of Bezier curve for each segment of the 
curve) 

ix. NGC_AII_Bezier_lnterpolated_values.m (Bezier 
interpolation of control points based on segmentation 
values) 

x. NGC_AII_Segments_Euclidian_distance_index.m (finds 
maximum square distance) 

xi. NGC_Bezier_lnterpolator.m (Bezier interpolation for given 
four control points) 

xii. NGC_Matrix_comparator.m (Compares and checks 
dimensions of the matrices) 

xiii. NGC_Bezier_CtrlPoint_getter (Least Square Method 
using specified Parameterization) 

xiv. plot_bezier_originaldata_controlP.m (plots curves and 
control points) 

xv. mdot_calc.m (calculates mass flow rate of the system) 


3. Outputs and generated curves from Matlab: 
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8m 


Pressure less 


Courbe de Bezier 



time [sec) 

Fig 11.4 

• End-user studies and evaluation of outputs: 

Pressure curve shows that the generated curves are converging, and 
mathematically it means the sequence has only one limit. In physics terms, it is 
attenuating (the gradual loss in intensity of flux through a medium). Converging 
results approves the validity of generated pressure curve from MATLAB 
program. 

B. Segment 2 (From manifold to roll block) 

1. Inputs: 

i. Characteristics of the components of the system (i.e. Manifold, 
dividing T, Diameter, Length, bend angle, surface roughness of 
the tube, etc.) 


Command Window 


^ To get started, select MATLAEi Help or Demos from the Help menu. 


enter 0 for Argon, enter 1 for Nitrogen: 0 
enter diameter of the tube in inch: .194 

e n t e r T emp e r at ur e at re gu 1 at o r in degree F : 

Fig 11.5 Command window in MatLab 
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ii. Mass flow rate, temperature, pressure and velocity of the flow 
at upstream calculated in segment 1. 

iii. NIST inputs 

2. MatLab files: 

i. mainR.m 

ii. DividingT_loss_coeff.m 

iii. Mach_interp_calc.m 

iv. Roll_loss_coeff.m 


3. Outputs: 

Mach2 

ratT2overT_0 

T_Roll 

ratC2overC_0 

c2 

rou_R 

U_roll 

P_roll 

Thrust_roll 

• End-user studies and evaluation of outputs: The analysis is finalized at this point and 
pressure and velocity are calculated at the roll block. For cross validation purposes the 
density calculated from the program, rou2=mdot_system/ (u_roll* A) should match 
the density from state of the gas provided by the NIST table: 



Fig 11.6 NIST program 

C. Segment 3 (From manifold to pitch block) 

Analysis is similar to segment 2, see section B for procedure. 
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Velocity (ft/s) Velocity (ft/s) 


Trends in Compressible adiabatic subsonic flow of gas 


Static Pressure 

Decreases 

Total Pressure 

Decreases 

Velocity 

Increases 

Density 

Decreases 

Temperature 

Decreases 

Mach number 

Increases 

Reynolds number 

Increases 

Stagnation temperature 

Constant 


Table 11.1 


Velocity in NMACS (Pitch) 



Pitch Velocity 


(52.6167 ft/s^p 


✓ 

/ 

■ 

* 

* 


/ 

* 

/ 


✓ 

_/ 

Tank Velocity 

— Regulator Velocity 

" (49.7 

h 

(50.134 ft/s) 


Station Number 
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i 
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V &♦*>*« Pimm 

V 
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Station Number 

Fig 11.7 

Synthesis of the Pneumatic System of Sounding Rocket 

In this section: 

• Break data into components of the pneumatic system by identifying causes and 
effects on overall efficiency of the system. 

• Synthesis and compile data together in a different way by arranging elements in 
a new pattern or proposing alternative solution. 


Pressure Loss in NMACS ( 


4000 


3000 


m Tank Pressure 
^\(3700 psi) 

\. 


2000 


1000 


Regulated Pressure Q* 
(700 psi) 


Roll Pressure" 
(177 psi) 

^ 


200 
190 
160 
140 - 
120 - 
100 
80 - 
60 


Velocity in NMACS (Roll) 


Roll Velocity • 
(194.34 ft/s) y 

✓ 

/ 

/ 

/ 

/ 

/ 

✓ 
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Tank Velocity 
J 1 (497 ft/s) 

tL 


_ (^Regulator Velocity 
m (50.134 ft/s) 


Station Number 
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Effect of key characteristics on pressure loss : 

Tube convergence, divergence, turns, surface roughness and other physical 
properties will affect the pressure drop. High flow velocities and / or high fluid viscosities 
result in a larger pressure drop across a section of pipe or a valve or elbow. Low 
velocity will result in lower or no pressure drop. Following table shows the effect of 
design factors on overall pressure drop of the system: 


factor 

Length 
increase 
by 1 
foot 

Diameter 
reduced 
to half 

90 

degree 

bend 

Surface 
roughness 
increased 
by 25% 

Manifold 

Preece 

valve 

Dividing 

tee 

Tube 

head 

loss 

0.25% 

more 

drop 

0.4% 

more 

drop 

0.10% 

0.05% 

4.50% 

73% 

3% 

0.38% 


Table 12.1 

Table 12.1 shows that maximum pressure drop happens in preece valves and in 
manifold. The pressure change is the sum of frictional losses and a pressure rise 
due to deceleration of the main pipe flow as some fluid is lost to the branch. 


Dividing tee and losses: 

The losses through a dividing tee are considerably reduced by rounding the 
edges of the T. in table below r is bend radius and D is diameter of tee. However the 
new manifold design eliminated the tee in the system. 


r/D 

0 

0.1 

0.2 

0.3 

0.5 

T Loss coeff. 

66% 

38% 

32% 

30% 

29% 


Table 12.2 


Tube surface roughness and material properties: 

Three different types of tubes have been analyzed: 

• Stainless Steel: 

a. Surface roughness is 40~70 micro inches. 

b. High yield strength preferable for vibration test and 
crack propagation criterion. 

c. Formability and bending issues. Stainless steel is 
hard to bend and the tube cross section issues. 
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d. Applicable for high pressure and temperature. 

• Aluminum: 

a. Surface roughness is 15 micro inches. 

b. Lower yield strength inappropriate for vibration test 
and crack propagation criterion. 

c. Aluminum is considered as semi brittle and easier to 
bend but it ripples. 

d. Aluminum handles pressure and temperature to 
certain degree. 

• Cupronickel (alloy 715 CuNi): 

a. Better surface roughness (10 micro inches). 

b. Copper Nickel alloy for high pressure and low 
temperature application. 

c. Better formability (copper tube, properly bent, will not 
collapse on the outside of the bend and will not buckle 
on the inside of the bend). 

d. Suitable for vibration test. 

Temperature of the operating fluid and behavior of the O-rings at specified temperature: 

Factors applying to all O-rings: 

• Compatibility 

• Temperature 

• Pressure 

• Extrusion 

Viton is a brand of synthetic rubber and fluoropolymer elastomer used in O-rings 
for pneumatic design. 

Redesign of the Manifold of the Pneumatic System (SolidWorks) 

The following design criteria are considered for the new manifold: 

• For achieving uniform flow the main inlet is sufficiently large in diameter 
and the main pipe flow velocity sufficiently small so that the pressure 
changes along the main branch are small compared with the pressure 
loss for fluid exiting through the branches. 

• Simple, robust design with standard mounting patterns for the valves, 
therefore; interchangeability becomes feasible for the system. 

• Coned shaped design at the end of the roll valve to reduce velocity 
therefore increase pressure. (Low velocity will result in lower or no 
pressure drop.) 

• Contains all the passages for an entire system with central shorter 
branches, and effective against external constraints. (Inlet and outlets are 
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located at appropriate surfaces for fewer bends in the tubes and it also 
eliminates the dividing T in the assembly.) 

Surface mounted compact polar design. 




Before 


After 



— 1.95 
— 
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Fig. 13.1 Pneumatic system Manifold 


Realization: 

Actual implementation of dynamic analysis details when implementing 
architecture isn't totally straightforward. Certain constraints are worthy of 
consideration: 

• Refining of discharge factor. C d is assumed as "most probable value" 
based on experimental data. However, it varies for different flow regimes 
(laminar, transient and turbulent). 

• System constraints: as mentioned in mathematical modeling section. 

• Locating of ACS system in the rocket and affect of the gyroscopic 
moves to interchange of energy. 

• Pneumatic component behavior: pneumatic components perform in a 
highly non-linear manner and the energy transfer medium is highly 
compressible. Both of these facts complicate the modeling and simulation 
of pneumatic system. 


Summary and Future Work 

Quite few analytical solutions exist for engineering problems. For instance in fluid 
mechanics, only simple potential flow has analytical solution. On the other hand, the 
real engineering problems usually are very complex and it’s impossible to solve them 
analytically and exactly. Hence, modeling and simulation is very powerful technique for 
real engineering problems. Simulation is used in this study in order to gain insight into 
functioning of the physical system. Key issues in this simulation include: 
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• Acquisition of valid source information such as measured temperature and 
pressure about the relevant selection of key characteristics and behaviors 
mentioned in this study, 

• The use of simplifying approximations and assumptions within the 
simulation, 

• Fidelity and validity of the simulation outcomes. 

For future trend, this project may conduct a further study to couple with 
application of the Kalman filter lor the system discretized in the time domain (time 
response modal analysis and Eigen value problem). 
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Appendix A 

Evaluation of derived equation by Buckingham n theorem: 




= R s .p.V~ 1 .A r . 

\ 



- 2 g(z r - z t ) 


In mathematical terms, if we have a physically meaningful equation such as 

f(x 1> x 2i ,...,x n ) = 0. 

The Xj are the n physical variables, and they are expressed in terms of 3 independent physical units, then 
the above equation can be restated as 

F{n 1 ,n2 , ... ,7 tp) = 0 

Where the Tij are dimensionless parameters constructed from the Xj by p = n - 3 equations of the 
form 


7Tj = X! m iX 2 m2 ... X- n mn 


Where the exponents m, are rational numbers (they can always be taken to be integers: just raise 
it to a power to clear denominators). 


*i = P(t) = ML -'T- 2 , x 2 =A 0 = L 2 ,x 3 = V - 1 = L~ 3 ,x 4 = RT c = - p = L 2 T~ 2 ,^ 

ML~ 1 T~ 3 ,x 6 = g = LT~ 2 , x 7 = (z r — z t ) = L 
F(7t 1 ,7t 2 , ...,n p ) = 0 -> 


Ap m 

&T 


->x 5 = 


% = ML -1 T -3 - ML~ 1 T~ 2 L 2 L~ 3 L 1 T~ 1 = ML- 1 ?- 3 - ML- 1 T~ 2 T-' l = ML -1 T -3 - ML - ^ -3 = 0 


The conclusion is derived Nonlinear differential equation is valid based on Buckingham n 
theorem. 
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Appendix B 

Gas Pressure recording and experimental data from pneumatic laboratory 

I. Non Flowing Gas recordings 

• Tank pressure - The initial pressure noted in the tank to perform a specific 
task, whether it be testing or the actual flight on the rail. The purpose of 
recording the tank pressure is to know what potential energy the system has 
in the form of gas to complete a specific task. This pressure is recorded by a 
transducer in the tanks manifold block that turns resistance into voltage 
ranging from 0-5 volts to give the pressure inside the tanks at any given 
time. 

• Manifold pressure - the supply of gas inside the manifold that incorporates 
the roll valves and feeds the pitch/yaw valves as well. The pressure here is 
important to the success of every mission and is directly controlled by the 
course regulator. To record this pressure there are two means of doing so, 
first is the transducer, similar to the transducer for the tanks reading and 
functions in the same manner, but for the final reading taking and for testing 
of the transducer. Readings are also taking from a direct gauge reading in 
which a gauge is screw into the main manifold to take reading during 
various phases in testing of pressure readings. 

II. Flowing Gas recordings 

• Tank pressure flowing is recorded by the transducer as in the non flowing 
readings and shows the supply of gas being used as well as the remaining 
gas. This is noted in most cases to check the efficiency, of the systems 
operation by the gas used at the end of the flight. 

• Manifold pressure flowing is recorded when gas is flowing to the valve or 
valves being activated. This reading is taking in two ways as well as the non 
flowing manifold pressure reading by the direct reading of the gauge 
screwed into the main manifold as well as the transducer reading. This 
reading is always recorded in the log of the system for the final flight 
parameters. 
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• Roll flowing pressure is taken by direct reading only, in where there is a tee 
pipe screwed directly into the roll block being fired and on the other end of 
the tee is the actual nozzle being used. To the side of the tee is the gauge to 
record the reading of the gas flowing. This reading is always recorded, and is 
a vital part of setting up the system flight pressures and nozzle 
configurations. 

• Pitch/Yaw pressure flowing - is the flow of gas from the pitch/yaw valve 
being fired, and is recorded in the same manner as the roll pressure. There is 
a tee installed in the block with nozzle on the other end and the gauge 
installed on the side of the tee. Along with the roll as well this plays a vital 
role in the flight pressure and nozzle installed into the system per flight 
requirements. 

. Lock up gas pressure 

• Manifold Lock up pressure is the pressure recorded whenever the valve is 
deactivated. This pressure is important for the valves and components in the 
system to insure the limits are within the proper specs. Similar to the 
manifold, and the manifold flowing pressure, this reading is taken directly 
from the manifold by the gauge being screwed into the block. This reading is 
also recorded and maintained in the logbook. 



Experimental dala from xducer reading 



Experimental data from xducer reading 



FigB.l 
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Appendix C 


Bezier Interpolation Technique 

Bezier curves were widely published in 1 962 by the French engineer Pierre Bezier, who 
used them to design automobile bodies. These curves were first developed in 1959 by Paul de 
Casteljau using de Casteljau’s algorithm, a numerically stable method to evaluate Bezier curves. 
Cubic Bezier curves defined by four control points PO, PI, P2 and P3 in either a plane or three- 
dimensional space were used for this study. The curve starts at PO going toward PI and arrives at 
P3 coming from the direction of P2. Usually, it will not pass through PI or P2; these points are 
only there to provide directional information. The distance between PO and PI determines ”how 
long” the curve moves into direction P2 before turning towards P3. 


P 


3 



FigC-1 

The first and last control points of Bezier curve are the first and last points of the input data 
segment. The input data can be separated into segments or just one segment by specifying the 
initial set of break points. But the second and third control points are determined by “least square 
method” . 

For fitting data, suppose we have a set of points (pressure drop in time domain extracted 
from Euler equation) and we want to approximate it using cubic Bezier. As an input we specify 
the value of limit of error (maximum allowed square distance between original and fitted data) 
and provide initial set of breakpoints. At least two breakpoints are required i.e., the first point 
and the last point of original data. Input data is divided into segments based on initial set of 
breakpoints. A segment is set of all points between two consecutive breakpoints. We have to fit 
each segment using cubic Bezier curve(s). Now the fitting process begins. 

We generate n points (approximated data) Q={q\,q2, . . . ,qn) using cubic Bezier 
interpolation such that cubic Bezier curve(s) passes through breakpoints. Then we measure the 
error between original and approximated (fitted) data. 

When approximated data is not close enough to original data i.e. limit of error bound is 
violated) then an existing segment is split (break) into two segments at a point called new 
breakpoint. After splitting the number of segments are increased by one (split segment is 
replaced by two new segments). Number of breakpoints are also increased by one (new 
breakpoint is added in the set of existing breakpoints). The point where the error is maximum 
between original and approximated data is selected as new breakpoint and this point is added in 
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the set of breakpoints. After splitting, repeat the same fitting procedure using updated set of 
segments and breakpoints until error is less than or equal to limit of error. 
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